BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, Float64[], tspan, Float64[])
Parsing parameters...done
Creating parameters...done
Parsing species...done
Creating species...done
Creating species and parameters for evaluating expressions...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
────
Time Allocations
─────────────────────── ────────────────────
────
Tot / % measured: 521s / 99.4% 208GiB / 99.8%
Section ncalls time %tot avg alloc %tot
avg
──────────────────────────────────────────────────────────────────────────
────
ODEProb SparseJac 1 456s 88.0% 456s 195GiB 94.2% 19
5GiB
ODEProb No Jac 1 32.0s 6.2% 32.0s 6.04GiB 2.9% 6.0
4GiB
Create ODESys 1 16.5s 3.2% 16.5s 3.76GiB 1.8% 3.7
6GiB
Parse Network 1 13.5s 2.6% 13.5s 2.13GiB 1.0% 2.1
3GiB
──────────────────────────────────────────────────────────────────────────
────
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show length(parameters(rn)) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 length(parameters(rn)) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = ModelingToolkit.varmap_to_vars(nothing, species(rn); defaults=ModelingToolkit.defaults(rn)) du = copy(u) p = ModelingToolkit.varmap_to_vars(nothing, parameters(rn); defaults=ModelingToolkit.defaults(rn)) @timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.) @timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.)
We also time the ODE rhs function with BenchmarkTools as it is more accurate given how fast evaluating f is:
@btime oprob.f($du,$u,$p,0.)
35.840 μs (3 allocations: 576 bytes)
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
────
Time Allocations
─────────────────────── ────────────────────
────
Tot / % measured: 1154s / 98.4% 217GiB / 99.7%
Section ncalls time %tot avg alloc %tot
avg
──────────────────────────────────────────────────────────────────────────
────
SparseJac Eval1 1 578s 50.9% 578s 5.35GiB 2.5% 5.3
5GiB
ODEProb SparseJac 1 456s 40.2% 456s 195GiB 90.2% 19
5GiB
ODE rhs Eval1 1 39.7s 3.5% 39.7s 4.01GiB 1.9% 4.0
1GiB
ODEProb No Jac 1 32.0s 2.8% 32.0s 6.04GiB 2.8% 6.0
4GiB
Create ODESys 1 16.5s 1.5% 16.5s 3.76GiB 1.7% 3.7
6GiB
Parse Network 1 13.5s 1.2% 13.5s 2.13GiB 1.0% 2.1
3GiB
SparseJac Eval2 1 114μs 0.0% 114μs 848B 0.0%
848B
ODE rhs Eval2 1 51.0μs 0.0% 51.0μs 752B 0.0%
752B
──────────────────────────────────────────────────────────────────────────
────
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)